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Abstract. We develop a numerical approach for computing the additive, multiplicative 
and compressive convolution operations from free probability theory. We utilize the reg- 
ularity properties of free convolution to identify (pairs of) 'admissible' measures whose 
convolution results in a so-called 'invcrtible measure' which is either a smoothly-decaying 
pg measure supported on the entire real line (such as the Gaussian) or square-root decaying 

^— I measure supported on a compact interval (such as the semi-circle). This class of mea- 

^^ sures is important because these measures along with their Cauchy transforms can be 

^N accurately represented via a Fourier or Chebyshev series expansion, respectively. Thus 

'H knowledge of the functional inverse of their Cauchy transform suffices for numerically 

^^ recovering the invertible measure via a non-standard yet well-behaved Vandermonde 

1*^ system of equations. We describe explicit algorithms for computing the inverse Cauchy 

QQ transform alluded to and recovering the associated measure with spectral accuracy. 

^' 
B 



o 



1. Introduction 



We propose a powerful method that allows us to numerically calculate the 'free' [21] ad- 
ditive, multiplicative and compressive convolution of a large class of probability measures. 
>- We see this method as complementing the symbolic techniques previously developed in 

[19] for so-called algebraic measures, i.e., measures whose Cauchy transform is algebraic. 



00 

in 

0\ Using the method developed in this paper, we can, for example, compute with spectral 

accuracy the free additive convolution of the semi-circle and the Gaussian which arises in 
[8] (see Figure [2]); or the free compression of the Gaussian which arises in [2] (see Figure 
CsJ |9|; or even the free additive convolution of the Gaussian with the counting measure on 

^ a single realization of a Gaussian Orthogonal Ensemble as a way to get insight on the 

^r*- rate of convergence to the asymptotic result in [8] (see Figures |2] & [3]). We go well 

S^ beyond these simple examples and hope that the proposed method allows practitioners to 

H experiment with these convolutions in the context of their experiments so they may find 

new applications of the underlying theory. 

We consider the free convolution operations on measures fiA and fiB (and compression 
factor a G (0, 1)) listed in the first column of Table [I] Each operation takes in one or 
two measures, and returns a new measure. What is known in each case is a relationship 
in transform space; i.e., there are transforms [211 [221 1211 E] R^iiy) and S^{y) so that the 
convolution operation can be represented simply as in the second column of Table [Tj 
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Operation Transform Operation Key Transform 



Free Addition 



G,.(.) = /'*'''"' 



Z — X 



xdM(x) 
FVee Multiplication s,,{y) = S.M^.M ^--' 



Free Compression „ 



v^ciy) = Rt^Ai<^y) 



G„(.) = /■*"<"' 



Z — X 

Table 1. Free convolution operations considered in this paper. 



The Cauchy transform of a measure /i on the real line is defined as: 

GJz) = / for z i. suppu. 

J z-x 

The key observation is that each transform i?^ and S^ can be expressed in terms of the 
functional inverse of the Cauchy transform G~^ (which we refer to as the inverse Cauchy 
transform). Therefore, we reduce the problem to the following two subtasks: 

(1) calculate the inverse Cauchy transform of the input measures pointwise; and 

(2) recover the output measure from knowledge of its inverse Cauchy transform. 

1.1. Invertible measures, their utility and the key underlying idea. In this paper, 
the class of admissible measures (see Section |3| are measures for which the inverse Cauchy 
transform (and hence the R or S transforms) can be accurately computed pointwise on 
an appropriate domain. 

The class of invertible measures, described next, are a subset of the class of admissible 
measures consisting of 

• square-root decaying measures: measures supported on an interval (a, b) with 
square root singularities at both the endpoints (such as the semi-circle), 

• smoothly decaying measures: smooth measures supported on the entire real line 
(such as the Gaussian), and 

• half square-root/smoothly decaying measures: supported on an unbounded inter- 
val (a, oo) or (—00,6) with a square root singularity at the finite endpoint. 

Invertible measures are a class of measures for which we can recover the output measure 
accurately from knowledge of its inverse Cauchy transform. The results of Section [2] 
state that the result of a free probability operation, if its support is simply connected, is 
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generically an invertible measure. The utility of the invertible measures can be discerned 
from Table [2l 

The key idea behind the proposed method is that invertible measures that are repre- 
sented via a Chebyshev or Fourier series expansion as in the second column of Table [2] 
have Cauchy transforms whose series expansions are closely related, as listed in the third 
column of Table [2j Thus, given a series truncation, we can efficiently compute the inverse 
Cauchy transform Gz^{yi) at yi by a companion matrix method. This knowledge of the 
inverse Cauchy transform G'^^iyi) at points {?/j}™^ coupled with the relationship (valid 
for y in the image of G^): 

G,{G-\y)) = y, 

implies that the desired series expansion coefficients {V'fcJiLi fo^ the Cauchy transform 
representation in the third column of Table |2] can be recovered by solving the Vandermonde 
system defined by: 

Gi,iG'^\y'd) = yi ioT i = l,...,m> n. 

This yields Algorithm [3] and Algorithm |5] for smoothly-decaying and square-root decaying 
measures, respectively. Choosing n appropriately yields the desired level of accuracy. 
Once these expansion coefficients are computed, we recover the measure fi via the series 
expansion in the second column of Table [2j 

The recognition that the class of invertible measures has a nice series representation 
for both the measure and its Cauchy transform is an important ingredient of the method; 
this insight originated in [T7] and might be of independent interest to free probabilists. 
Representing the measures via another basis that yields a sparser series representation of 
the measure but that does not yield a sparse, directly computable and invertible, series 
representation of Cauchy transform does not lead to an algorithm for computing the 
inverse Cauchy transform, thereby stalling progress. 

The paper is organized as follows. In Section [2| we discuss the analytic properties 
of the Cauchy transform that we use to construct the numerical method. In Section |3| 
we describe a numerical approach for the first sub-task, i.e., the calculation the inverse 
Cauchy transform for several types of admissible measures that arise in practice. We 
then solve the inverse problem in Section |4j we develop an algorithm to recover an 
unknown measure based on pointwise evaluation of its inverse Cauchy transform. In each 
stage, we achieve spectral accuracy, whenever we know the form of the measure. In the 
remaining sections, we apply this numerical algorithm to free additive, multiplicative and 
compressive convolution. 

2. Regularity properties of free convolution and its implication 

We think of admissible measures as candidate 'input' measures that we would like 
convolve using the operations in Table [Tj In this viewpoint, invertible measures are the 
generic 'output' measures that result from the convolution of admissible 'input' measures. 

Table |2] lists the class of invertible measures; recall that these are measures that can 
be recovered accurately from knowledge of their inverse Cauchy transform. In contrast, 
admissible measures are those for which we can compute the inverse Cauchy transform 
accurately. 
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Type 



Measure 



Cauchy transform 



Square- Root 

Decaying 

(e.g. Semi-Circle) 



afj,(x) = ip{x) : ax, 



b — a 



i^i^) = T.T=o^kU^{M^^\^{x)) 



where J , (z) = z — \J z^ — 1 and, 



A^(a,6)(2) ■ 



2 "•" 2 ^ 



dii{x) = ip{x) dx, 

1 + x — a 



Half Square Root 
/Smoothly Decaying '^^'^ = ^^=o^>^U,(M-]^^{x)) 



where J7 (z) = z — ^/z^ — 1 and, 



A^(a,oo)(^') =a+ T^ 



l + x 



Smoothly 

Decaying 

(e.g. Gaussian) 



dii(x) = i/>(x)dx, 

where ij})^ = ^_j. and, 
M(a;) = ^ 






-27r • 



Efe^-l^fc«w^3w<o ' 



where u{z) ■ 



Table 2. Series representation of invertible measures and their associated 
Cauchy transforms. 



The semi-circle and Gaussian measures are both invertible and admissible; the uniform 
measure on an interval and the (discrete) point measure are admissible but not invertible. 
Invertible measures are thus a proper subset of the class of admissible measures. Might 
this be a shortcoming of our proposed method? We assert otherwise and argue why the 
mathematics of free convolution gives us license to carve out the smaller class of invertible 
measures from the larger class of admissible measures. 

Simply put, the free convolution of two admissible measures generically results in an 
invertible measure. An important implication is that we can predict the form of the con- 
volved measure and apply a suitable algorithm (see Section H) for recovering the measure 
from its inverse Cauchy transform. 

In this paper, we restrict ourselves to measures that have a density; let /xi be strictly 
positive everywhere then for any non-trivial ^2 we have that /xi ffl ^2 is absolutely contin- 
uous with an analytic density supported on the entire real line. This follows directly from 
the results in [3]. The subordination theory of free convolution implies that the worst 
behavior of the two measures is preserved in the convolution. These statements carry 
through for free multiplicative convolution as well. Thus the convolution of admissible, 
smoothly decaying measures results in an invertible smoothly decaying measure. 

Now consider two (non-point) measures; each supported on an interval (that decay 'fast 



enough' at the edge of the support). Then, by Theorem A.l, the resulting measure will 



be an invertible square-root decaying measure (thanks to Serban Belinschi for supplying 
the proof) and supported on an interval (see |2l])- A similar argument (omitted here) 
holds for free multiplicative convolution as well. 
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The half-square root /smoothly decaying invertible measure arises from the free convo- 
lution of a measure supported on an interval with a measure supported on a semi-infinite 



interval. The argument in Theorem A.l may be adapted to show square-root decay at 
one edge of the support. 

The setting that is the least well-characterized in this framework is the convolution of 
two point measures. Here anything can happen. For example, the convolution of /ii(x) = 
0.5S{x) + 0.56{x — 1) with fi2{x) = /ii(a;) yields the arc-sine distribution which is neither 
admissible nor invertible (according to our definition). An improved characterization of 
the class of measures that results from the free convolution of point measures would help 
determine if and how the proposed method may be extended to this setting. 



3. Computation of inverse Cauchy transforms 

We consider the numerical calculation of the Cauchy transform and its inverse for four 
types of measures, which we refer to as admissible: 

(1) smoothly decaying measures: smooth measures supported on the real line; 

(2) square root decaying measures: measures supported on an interval (a, b) with 
square root singularities at the endpoints; 

(3) half square root/smooth decaying measures: measures supported on (a, oo) or 
{—oo,b) with a square root singularity at the finite endpoint; 

(4) point measures and 

(5) combinations of the above. 

For brevity, we omit the details for half square root /smoothly decaying measures below, 
as they can be treated very similarly to square root decaying measures. We include the 
relevant formulae in Table |2l 

We note that the formulae for the Cauchy transforms below follow from Plemelj's lemma 
[13] : i.e., if d/i = ipdx for suitably smooth ip, then 

(j)'^{x) — (p^ix) = —2'Kiil){x) and 0(oo) = 

if and only ii cj) = G^, where 0+ denotes the limit in the complex plane from above and 
(i>~ denotes the limit from below. 



3.0.1. Computing the Cauchy transform and its function inverse of smoothly decaying 
measures. A smoothly decaying measure is of the form 

d/i(x) = ip{x)dx, 

where tp G C^(— oo, oo), ip' has bounded variation and tp{x) = - + 0(x~^) as a; — )■ ±oo 
for some constant a. We can expand 

oo y . X A: 

k=—oo ^ ' 
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where ipk = '4'-k (since ip is real-valued) and '?/'(oo) = J2 
transform satisfies [271 E! 



oo 
fc=— oo 



(-)V. 



0. The Cauchy 



G,iz) 



-2tii 



2-^k=—oo ^k \ij^z) 



Q^Z > v^ , .u , 



(2) 



If?/' is C°°(— cxD, cxd) and ip has a full asymptotic expansion that matches at ±cxd, then the 
series ([l]) converges spectrally quickly. Moreover, we can rapidly compute the coefficients 



of the expansion by applying the FFT to the pointwise function samples ip ( i 1 1 u'" ) ? where 
Mm are m evenly spaced points on the unit circle: 



Ur 



■l,e^ 






Thus we take m = 2n + 1 and approximate 



G,{z) 



-27Ti 



For large n, T/^fc are accurate to machine precision. 

Now consider the problem of computing Gt^^. Note that 

l-z 



k=0 



'i^k 



G,. 



l + z 



-27ri 



ELo^kz' \z\<0 

-E'kLn^kz' \z\>0 



n 

E( 

fc=0 



'^k 



We can therefore invert the approximation of G^ using a companion matrix method. In 
detail, we compute the eigenvalues {Xf (y) , . . . , A.^(y)} of the matrix 



/ 



V'0-E^^0(-)''^fc-3|^\ 



V 



1 



■4>n 
■01 



V'n 



/ 



Similarly, we compute the eigenvalues {A^^ (y), . . . , A„ (y)} of the matrix 



/ 



i>-n 



\ 



Then 



where 



G-'iy) 






^ 1 + A(y)' 



A(y) = {A+(y) : |A+(y)l < 1} U {A7(y) : |Ar(^)| > 1}. 
Because the Cauchy transform is single valued at infinity, there is a neighbourhood of 
zero such that, for large n, X{y) is single valued. 
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3.0.2. Computing the Cauchy transform and its function inverse of square root decaying 
measures. Suppose that n is supported on the interval (a, b) with the form 

2y/x - ay/b - X 

djjix) = tpix) da;, 

b — a 

where ip & C^ [a, b] and ip' has bounded variation. We can represent 

oo 
fc=0 

where Uk denote the Chebyshev polynomials of the second kind and M(^a,b) is an afiine 
transformation from the unit interval to (—1,1): 

a + b b — a 
Then the Cauchy transform satisfies 

oo 

G,(^) = 7^5^^,_lJ;l(Af(;^,)(^))^ (3) 

fc=i 
where 

J-l(2;) = z- y/z- iVl + Z 
is an inverse to the Joukowsky transform 

Remark 1. We have not found this exact expression for the Cauchy transform of a square 
root decaying measure in the literature, though directly related expressions are in pTllTS] . 
In short, it follows from Plemelj's lemma and the fact that 

hm i^ — = (7fc_i(x)Vl - x^, 

verifiable by the substitution x = cos 9 [16] . 

We can compute the coefficients ipk whenever we can evaluate ip pointwise, and thence 
the Cauchy transform itself. This is accomplished by first computing the expansion in 
terms of Chebyshev polynomials of the first kind 

n-l 

V'(M(„,fe)(x))^^^fcTfc(x), (4) 

fc=0 

which can be accomplished by applying the DCT to '?/'(M((j f,)(x„)), where x„ are n Cheby- 
shev points of the second kind: 

We then transform the expansion (H) to a expansion in terms of Chebyshev polynomials 
of the second kind using the formulae 

To(x) = Uo{x),T,ix) = ^ and T,(x) = ^'^""^ - Uk-2{x) ^^^ ^ ^ 3, 3, . . . 

This approximation will converge spectrally when ip G C°°[a,6]. 



8 SHEEAN OLVER AND RAJ RAO NADAKUDITI 

3.0.3. Computing the inverse Cauchy transform of a square root decaying measure. We 
want to solve 

G^,{z) = y. 

Since J^^{J{w)) = w ior w inside the unit circle, we have 

n n 

k=l k=l 

We can thus solve G^ (M(a^fe)(J(w))) = y to find w{y) inside the unit circle using a 
companion matrix method (as above). Then 

G-^\y) = M^aMJiwiy))). 

3.0.4. Computing the Cauchy transform and its function inverse of a point measure. Sup- 
pose d/i(x) = 6{x — a)dx. Then its Cauchy transform is trivial: 

J z — X z — a 
Its inverse is 

3.0.5. Combination of previous types of measures. Consider the case where /i is a sum of 
point measures, for example, the counting measure 



1 " 
du = — > 5(x — \i)dx 



n . 

of one realization of a n x n random matrix with eigenvalues Ai, . . . , A„. The Cauchy 
transform can be computed directly using the previous approach, however, its inverse is 
no longer straightforward to compute. To calculate the inverse, we surround the support 
of fi by an ellipse E(^a,b),r in the complex plane, on which the Cauchy transform of the 
measure is smooth. We then exploit analyticity of the Cauchy transform outside of this 
ellipse. 

Define an ellipse E(^a,b),r surrounding the interval (a, b) as the image of the unit cir- 
cle under the map M(a,b)iJii"w)), with inverse Kj^^ (Mr^\.(z)) . We can then expand a 
function g defined on Ei^a,b),r by 

oo 

g{M^aMJi^^))) = E 9kw', (5) 

fc=— oo 

where the coefficients are computable numerically using the FFT as before. 

On and outside this ellipse, G^ is analytic and vanishes at infinity, therefore G'^(M((j {,)( J(rw)) 
is analytic inside the unit circle for r < 1 and vanishes at zero. Hence we can efficiently 
represent it in terms of positive Laurent series: 

oo p-. -| fc 

fc=i L 
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Analyticity of this sum implies that the expression holds true for z outside -E{a,fe),r as well. 
Mapping this sum back to the union circle allows us to compute G~^ using companion 
matrix methods. 

Remark 2. Note that r is a free parameter. As r approaches one, the ellipse approaches 
the support of the measure. Since G^ generically has singularities on the support of /i, the 
convergence rate of the expansion ([s]) degenerates. For r small, the ellipse is too large and 
the region of validity for computing G~^ shrinks. For the numerical examples below, we 
fix r arbitrarily (r = .8). A better approach would be to exploit the connection with the 
closely related problem of optimizing the radius of circle used in numerical differentiation; 
a problem solved in [^ . 



4. Recovering a measure from its inverse Cauchy transform 

Using the preceding formulae and the expressions for the transforms below, we can 
successfully compute the inverse Cauchy transform Gt} of some unknown measure [i 
pointwise, which will arise as the output of a free probability operation. If /i is either a 
smoothly or a square root decaying measure, we assert that, under broad condtions, the 
following algorithm will construct an accurate approximation to /x: 

Algorithm 1. Compute measure from inverse Cauchy transform 

Given G^{y) which is defined pointwise for y E G^{C), point cloud jm = {Vi, ■ ■ ■ , Vm) 
in the upper half plane and the assumed form of the measure /i (smoothly or square root 
decaying measure); compute a representation /i as follows: 



1: Use Algorithm p^ to prune y^/ so that all points lie inside G^{£)] 

2: If the desired form for /i is a smoothly decaying measure, use Algorithm |3| 

3: Otherwise, if the desired form for /i is square root decaying measure, use Algo- 
rithm m 

The first step of the algorithm is to assure that all sample points lie within G'^(C). 
Algorithm 2. Prune points 

Given G~^{y) which is defined pointwise for y E G^{C) and point cloud ym = {vi, ■ ■ ■ , Vm)', 
compute Ym C G'^(C) as follows: 

1: Select the elements of y*/ that satisfy sgnQy ^ sgnQG'^^{y). 

Proposition 4.1. Suppose that fi is smoothly or square root decaying measure. IfGz^{y) 
can he analytically continued in a domain containing G^{C), then there is a domain D 
containing G^iC) such that for all y E D, ^y j^ 0, sgnQG'^^{y) ^ sgnSy if and only if 

yeG.iC). 
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Proof. We first show that G~^ has no turning points (points where (G^^)' vanishes) in 
G^i(C) except possibly at Gfj,{a) and G^{b). Note that 

From (|3| the only points where G' can blow up are at the endpoints (a, b) for square root 
decaying measures. Similarly, from (pi) the only point where G' can blow up is at infinity 
for smoothly decaying measure. 

Because of the positivity of measures, we know that '^Gfj_{z) is negative for 5^2 > and 
positive for '^z < 0. Therefore, for y G (j'^(C), sgnQG~^{y) ^ sgnQ^y. Moreover, since 
G^^ is analytic in a domain containing G^(C), and since no point in the complex plane 
on G^(supp/i) is a turning point, it is true that sgnQ=G^^(y) = sgnQ^y ior y E D outside 

G,(C). 

D 

If we assume the measure is smoothly decaying, then we know precisely the form of its 
Cauchy transform, but we do not know the relevant coefficients of the expansion. The 
following algorithm computes these coefficients by applying least squares to the equation 

G,{G-^\y)) = y, 

which is valid for y G G'^(C). 

Algorithm 3. Compute smoothly decaying measure 

Given G~^, point cloud y^ inside Gf^{C)r\{z : '^z > 0} and positive integer n; compute 
a representation of fi that is smoothly decaying as follows: 

1: Compute ipk by solving the following system in a least squares sense: 



-27vi 



£-(^^)'-&->^- 



Vj- 



Define i/jq = -2^Y.l=ii~^)^i'k and tp-k = ipk- Then 

d/i ^ V" V^fc ( — — I dx. 



Theorem 4.2. Suppose that fi is a smoothly decaying measure, and {ym} o-i^^ o- sequence 
of sets of m points lying inside G^{C) H {z : '^z > 0} which cover G^(C) H {z : Qz > 0} 
as m ^ oo at a sufficiently fast rate (see proof and Appendix B for precise definition). 
Then there exists m sufficiently large depending on n so that the Algorithm [^ converges 
to /i as n — 7- oo. 

Proof. Because of symmetry, including y™, in the least squares system will not alter the 
approximation of /i. Therefore, denote [ym,ym] = [vi, ■ ■ ■ ,y2m]- Then 
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for some (unknown) Zj inside the unit circle. Under this transformation, the least squares 
system takes the form 

n 

fc=i 
This is a Vandermonde system, with an unusual distribution of points. However, as 
m — )■ cxD, the points Zj must cover the unit circle, and therefore convergence follows from 
Lemma IB.2I 

D 

Remark 3. Though we omit the details, it is straightforward to prove that the point sets 
Ym that we generate below satisfy the conditions of the preceding theorem, due to the 
analyticity of the operations involved. 

We can adapt this approach to square root decaying measures as well; since, assum- 
ing we know the support of the measure, we again know a precise form for its Cauchy 
transform. 

Algorithm 4. Compute square root decaying measure 

Given G^^, point cloud y^ inside G^(C)n{z : 3z > 0} and positive integer n; compute 
a representation of /i that is square root decaying as follows: 

1: Compute (a, 6) ~ supp/x using Algorithm [SJ 

2: Compute (real- valued) ipk by solving the following system in a least squares 

sense: 

-nY^i^u-i^J-' (^K6)(G;'(y,)))' ^ %• and 
fc=i 

k 



where 



vr^^,_i3J;MM(;^,)(G,i(y,)) ^%, 



fc=l 



a + h b — a 



Then 



k=0 



, Zy/X — ay/0 — X JT-^ , ^. , . . , -., 

d/i ^ J— ^ipkUk{M(a,b){x))dx. 



If supp /i is calculated accurately, the convergence of Algorithm |4] follows by the same 
logic as Theorem |4. 2 Thus we are left with one last task: computing supp/i. 



Algorithm 5. Compute the support of a square root decaying measure 

Given G~^, its first two derivatives and initial guesses (ao,&o); compute an interval 
(a, b) approximating the support of fi as follows: 
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Figure 1. A plot depicting the curves on which G J^ is real. These include 
Ff, the image of G^, and F^, the image of G~ and (^a^^fe), for ^a = G^{a) 
and ib = G^{h). 

1: Compute a and h by solving G~^ {y) = using Newton iteration, with oq and Bq 
as initial guesses. 

While we only discussed the computation of G~^, computing its derivative is straight- 
forward since 

^^"'^'^^^ ^ G',{G-^^{y)) 
and the derived formulae for Cauchy transforms of admissible measures can be trivially 
differentiated. Similar logic allows us to compute the second derivative needed to perform 
the Newton iteration. 

The convergence of the algorithm follows: 

Proposition 4.3. Suppose fi is a square root decaying measure. For a sufficiently accurate 
initial guess, Algorithm\^will converge to supp/i. 



Proof. Let [a,b] = supp fj, and 



F± = G^([a,6]) 



on each of which Gz^{y) is real valued (as its image is [a, h\). The Cauchy transform is real 
valued at the endpoints of the support (since the integrand of the Cauchy transform is 
integrable), therefore V± intersect the real axis at ^a = G^{a) and ^b = G^{b), as depicted 
in Figure ni Moreover, G~^{y) itself is real valued on the real axis (since the Cauchy 
transform is real valued on (6, oo) and (— oo,a)). Therefore, ^a and ^b are saddle points 
of Gz^{y), and the Newton iteration is guaranteed to converge for sufficiently accurate 
guesses. 

D 

5. Free additive convolution and the R transform 

The R-transform, defined as 

R,iy)--=G,\y)-l/y, 

is the analogue of the logarithm of the Fourier transform for free additive convolution. 
The free additive convolution of probability measures on the real line is denoted by the 
symbol ffl and can be characterized as follows. 



FREE PROBABILITY CALCULATOR 13 

Let An and B^ be independent n x n symmetric (or Hermitian) random matrices that 
are invariant, in law, by conjugation by any orthogonal (or unitary) matrix. Suppose 
that, as n — > oo, fiAn — ^ f^A and fiB„ — > fJ'B- Then, free probability theory [211 El] 
states that fiA„+B„ — > At^ H /^b, a probability measure which can be characterized in 
terms of the i?-transform as 

^MAfflMs (y) = ^MA (y) + ^Mfl (y) ■ (6) 

Rearranging ([6]), we find that 

G,l^,,{y) = G-^liy) + G;^l{y) - I 

Therefore, if ^a and fiB are known in admissible form, we can apply the methods of 
Section p^ to compute G~^^^^. To employ the algorithm of Section E, we need to assume 
the form of /i^i ffl yU^- The discussion in Section |2] allows us to know this apriori. 

Our last task is to generate a point set ya/ that covers (^^^^^^(C). To accomplish this, 
we use the following theorem: 

Theorem 5.1. 

G,,ffl,,(C)cG,,(C)nG,,(C). 

Proof. This statement first appeared in Proposition 4.3 of [26]. It is a direct consequence 
of the subordination [6] of the functions expressed in ([T]). D 

Thus we generate jm by first generating a point set z^^^m off the support of fiA and 
then use the point cloud ym = G'^^(z^^^m)- Because we impose that the measures are 
real, we automatically have that G^{z) = Gfj,{z). Therefore, as required in the algorithms, 
we restrict z^^m to the upper half plane. Many suitable approaches for generating z^^m 
exist, the approach we take is the following: 

Algorithm 6. Generate point sets 

Given an a smoothly or square root decaying measure /x; compute a set of points z^^m 
lying off supp /i as follows: 

1: Generate a point cIa/ on the unit disk by taking a tensor product of u^ with 

Mz\Jxm), the m Chebyshev points on (0, 1); 



If supp/i is the real line, return ij^ 



l-dM . 

M ' 



3: If supp/i is an interval (a, b), return the subset of M(^a,b){J{dM)) which lie in the 
upper half plane. 



5.1. Numerical examples. 

Remark 4. Throughout th( 

distributions unless otherwise specified. Sn denotes an n x n random symmetric matrix. 



Remark 4. Throughout the paper, we use mean zero and variance 4= for Gaussian 
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Figure 2. Free addition of a Gaussian distribution with a semicircle distribution. 
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Figure 3. Free addition of a Gaussian distribution with a single instance 
of an approximate semicircle distribution. 



constructed by generating a random matrix An with Gaussian distributed entries and 
defining 



^n 



An + A^^ 

\p2Jn 



Qn denotes a random orthogonal matrix, generated by computing the QR decomposition 
of An- Finally, we generate a histogram associated with a random matrix ensemble i?„ 
by computing the eigenvalues of 100 instances of i?„. 

In Figure [2| we plot the numerically calculated free addition /i^ ffl /U5 of a Gaussian 
distribution /ic with a semicircle distribution ^s- This distribution was shown in |S] to be 
the limiting eigenvalue distribution of a class of Markov matrices. The left graph contains 
a plot of a Gaussian distribution (dotted), semicircle distribution (dashed) and their free 
addition (plain). The right graph compares the computed free addition with a histogram 
of Qi5oAi5oQ75o + 'S'150, where A„ is a n x n diagonal matrix whose entries are Gaussian 
distributed. 

Often one does not have exact expressions for the limiting distributions of the eigenval- 
ues, but rather, one can sample a single instance from the distribution. In this case, the 
counting measure — a sum of point measures — over this single instance can be calcu- 
lated. In Figure [3] we repeat the experiment of Figure [2] where the semicircle distribution 
is replaced with the counting measure /iAgo of a single matrix ^50 drawn from 6*50. On 
the right, we compare the computed distribution with the histogram of QsoA-soQjo + ^50) 
where ^50 is now a fixed matrix. 
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Figure 4. Free addition of a Gaussian distribution with approximate semi- 
circle distributions for n = 100, 200, ... ,400 (left). The scaled (by n) 
Kolmogorov-Smirnov distance (D„ = sup x\Fn{x) — F{x)\ where F{x) is 
the distribution in Figure |2] and -F„ is the distribution in Isl) between the 
cdfs illustrating convergence in the respective cumulative distribution func- 
tions (right). 
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Figure 5. Free addition of a Semicircle distribution with the equilibrium 
measure associated with the potential V{x) 
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As n — )■ oo, fis^ ffl fiG will converge in some sense to /^^ ffl /ig, as seen in the right hand 
of Figure |4j We can estimate this growth by comparing the maximum difference of the 
cdf of computed measures for growing values of n. In the right-hand side of Figure 4| we 
plot this scaled by n, demonstrating that the convergence rate appears to be 0{n~^). 

In Figure |5] we compute a measure which is square root decaying. Here we define 
/i4 as the equilibrium measure of the potential V{x) = x^ (see [20] for definition of 
equilibrium measures), which we know in closed form [TU]. We then calculate /i5ffl/i4 using 
Algorithm [1} There is no obvious way of generating a histogram for this measure; hence, 
unlike other examples, there does not exist a Monte Carlo approach for approximating 
/i5 ffl /i4. However, this is an example which was calculated symbolically in |19j, hence 
we can compare our numerically computed measure with the exact measure. We plot the 
error for ra = 20 (dotted), 40 (dashed), 60 (dash-dotted) and 80 (plain) as M increases. 
Recall that n is the number of coefficients in the Chebyshev representation of ^s and 
/i4 while M is the number of points in the point cloud used in the least-squares based 
measure recovery algorithm described in Algorithm 6. The error is computed by taking 
the maximum error over 100 Chebyshev points on the interval supp(/is ffl yU4). 
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Figure 6. Free addition of a Semicircle distribution with the equihbrium 
measure associated with the potential V^(x) = e^ — x. 
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Figure 7. Free addition of a step distribution with a semicircle distribution. 



In Figure |6| we define \iem as the equilibrium measure of the potential V{x) = e^ — x 
— which we calculate numerically (in the required form) using the approach of [18] — 
and then calculate i^s^I^'EM- This is an example which cannot be computed symbolically, 

as in [T9]. 

Finally, in Figure [7] we calculate the free addition of a semicircle distribution with a step 
distribution /i^ ffl (2l(-i,i)), by assuming that it is square root decaying, demonstrating 
that this behaviour is, in fact, generic. While |l(_i.i) does not have the form considered 
in Section [3} we can in fact calculate its Cauchy transform and inverse Cauchy transform 
explicitly: 



2l(-l,l)^ ' O 



G\ 



-1 



{y) 



coth I + tanh ^ 



Therefore, the algorithms of Section |4] are still usable. We compare the computed distri- 
bution with the histogram of Qsoo^sooQjoo + 'S'soO) where A„ is a diagonal matrix whose 
entries are evenly distribution on (—1, 1). 
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6. Free multiplicative convolution and the S transform 

In the case where fi ^ Sq and the support of /i is contained in [0, +oo), one also defines 
its T -transform 

/x 
d/i(x) for z 4: supp/i. 
z-x 

The S -transform, defined as 

S,{y):={l + y)/{yT-\y)), 

is the analogue of the Fourier transform for free multiplicative convolution Kl. The free 
multiplicative convolution of two probability measures fiA and fiB is denoted by the sym- 
bols Kl and can be characterized as follows. 

Let An and Bn be independent nxn symmetric (or Hermitian) positive-definite random 
matrices that are invariant, in law, by conjugation by any orthogonal (or unitary) matrix. 
Suppose that, as n — )■ oo, fi^,, — >■ jj^a and /iB„ — > /^b- Then, free probability theory 
states [22] that /iA„-B„ — > fJ'A^ fJ'B, a probability measure which can be characterized in 
terms of the S'-transform as 

Sfx^mnBi^) = Sf_i^{z)S^g{z). 

The T transform can be computed in the same way as the Cauchy transform; we only 
need to multiply the function ip by x beforehand. This is true of its inverse as well. From 
the relationship of the S transform, we know that 

Note that T^^^/j.^ = G^^, for the (non-probability) measure fie defined by 

dnc{x) = xd[fiA^lJ'B]ix). 

Therefore, we can use the Algorithm [l] to find d/ic, and in turn fj,A Kl fj,B. Similar to free 
addition, we use the point cloud yM = ^m^(z^^,m)- 

Theorem 6.1. 

T,,^,,{C)cT,,{C)nT,,{C). 

Proof. This is a direct consequence of the subordination of the functions expressed via 
the relationship T^^(F4(z)) = T^^^^^(2;) = T^g^Fsiz)) derived in [HI Theorem 3.5, pp. 

157]. D 

In Figure [8} we consider the problem of computing a free product of a a shifted semicircle 
distribution with a singular Marcenko-Pastur distribution 



dfiMp{x) = — — ^dx. 

While this distribution is not admissible, it is when we multiply by x; as in the definition 
of the T-transform. The procedure then works as before. We compare the computed 
measure with a histogram of 

-B200-B2oo('S'200 + 3/), 

where Bn = T^^n and An is a.n n x n random matrix with Gaussian distributed entries, 
now with mean zero and variance one. 
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Figure 8. Free times of a shifted semicircle distribution with a Marcenko- 
Pastur distribution. 

7. Free compression 

Let Bn be the n x n matrix generated by the taking the upper left n x n block of 
Qm-^mQlni where n < m. If the eigenvalues of Am tend to the distribution /x, then the 
eigenvalues of Bn tend to the free compression of /i, i.e., 

m 
n 

Let a e (0, 1]. We have that [H] 

RaB^l{z) = R^{oiz). 

Rearranging the definition of the R transform, we find that 

Therefore, we can apply Algorithm [T] to compute a H /x, with the point cloud yu = 



7.1. Numerical examples. In Figure^ motivated by the theoretical results in [2], we 
compare the compute free compression of a Gaussian distribution with a histogram of the 
a300 X a300 principal block of Qsoo-^sooQsoo, where A„ is an n x n diagonal matrix whose 
entries are Gaussian distributed. 



8. Extensions 

We now identify some extensions of the proposed method: 

• Measures supported on multiple intervals. Here there are two issues that must 
be overcome: computation of the inverse Cauchy transform, and determination of 
the support of the measure. The major complication is that the inverse Cauchy 
transform is multi-valued. 

• Free rectangular convolution (see [5]). This operation inherently requires compu- 
tation with measures supported on multiple intervals. 

We plan to release a software implementation based on the ideas described in this paper. 
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Figure 9. Four compressions of a Gaussian. 
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Appendix A. Square-root decay as generic edge behavior of a convolved 

MEASURE 

Theorem A.l. Assume that both fii and fi2 ore compactly supported. If fii and ^2 do 

not have 'very fast decay' at the boundary of their respective supports, then fii ffl /i2 is 
square-root decaying at the left and right edge of the support. 

Proof. We restrict ourselves here to probability measures /xi and /i2 which do not have 
'very fast decay' at the boundary of their respective supports. The meaning of 'very fast 
decay' will be made a bit more precise, but not completely precise, in what follows. Let 
yU3 = /ii ffl /i2 and for j = 1, 2, 3, let Fj{y) = l/G^^ (y) denote the reciprocal of the Cauchy 
transform of the measure fij and let 0j denote the subordination function: 



F,{0,{z)) = F2{02{z)) = Fs{z), 



(7) 
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where by the resuhs of Bercovici and Behnschi pQ, we have that the function 02 (note 
01 (z) = F2 (02(2:)) — 02 (^) + z) satisfies the fixed point equation: 

02{z) = Fi (F2 (02(2;)) - 02(2) + z)- (F2 (02(2;)) - 02{Z) + z) + z. (8) 



The right hand side of ([8]) can be expressed as a function of two variables f{z, w) = 
Fi{F2{w) — w + z) — {F2{w) —w + z) + z evaluated at w = 02(2;). It exhibits good behavior 
and maps the product of upper half-planes to the upper half-plane. However, what is 
interesting for our purpose is the fact that 02 fails to be analytic at a point a;o € M if and 
only if: 

• F[ (F2 (02(xo)) - 02(xo) + xo) {F^ (02(xo)) " 1) " i^2 (02(a;o)) = or, 

• one of Fj fails to be analytic at 0j{xo). 

This second possibility, as can be seen from the first condition (recall that 0i{z) = 
-^2(02(2;)) — 02{z) + z), requires that the product {F[{0i{x)) — l)(F2(02(a:)) — 1) stays 
strictly less than one when 0j{x) moves all the way to the support of the continuous part 
of fij] this is the sense in which 'the decay at the boundary should not be very fast'. 

Now assume that this is the case and that {F[ (0i(x)) — 1) (Fg (02(x)) — 1) = 1 has a 
solution at a point xq G M for which both 0j{xo) are real and in the domain of analyticity of 
Fj. Thus 02 is analytic and we are in a position to apply Weierstrass' preparation theorem 
(see, for example P]) to the function {z,w) H- f{z,w) — w around {z,w) = (xq, 02(a^o))- 
The theorem states that: 

f{z, w)-w = (w^ + ci{z)w^-^ + ■■■ + Ckiz))(j){z, w), 

where (f){z, w) is nonzero on a neighborhood of the chosen point, Cj are just analytic maps 
around Xq, and k is the number of zeros with multiplicity of w H- f{xo,w) — w on some 
small enough neighborhood. 

An application of this theorem reveals that right of the supports of the /ij's, we have 
that F'-{w) < 0, F'jiw) > 1, and we obtain that the second derivative of w H- /(xq, w)—w 
is strictly negative at w = 02(3^0)- 

Thus, k = 2, and 02 has a square root singularity at xq. Applying the same argument 
we can also show that 0i has a square root singularity at Xq, and hence so does F3 
since F^i^z) = Fi{0i{z)) = -^2(02(2;)). This implies that G^^{z) will have a square root 
singularity at Xq, the endpoint of the support and hence, by Plemelj's inversion formula, 
/i3 will be square-root decaying. The same argument works for the left-hand side of the 
support as well. D 



Appendix B. Convergence of Vandermonde systems with large number of 

POINTS 

The following proofs are straightforward (we thank Ben Adcock for help proving them), 
though we have not found them in precisely this form in the literature. 
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Proposition B.l. Let d^ = {di, . . . ,dm) be a point cloud that covers the unit disk as 
m — )■ oo, and 



(\ di ■■■ rf^-i^ 



V 



1 d„ 






the mxn Vandermonde matrix associated with the point cloud. Then for any < e and n, 
there exists m large enough so that \\V^\\ < ^/n + 5, where V'^ denotes the Moore-Penrose 
pseudoinverse of V . Furthermore, if, for all \x\ < 1 and e > 0, the smallest m such that 
min(|a; — dm|) < e satisfies — = 0(e") for some a > 0, then m = 0{n°'). 

Proof. Assuming that V^ has full column rank (which will follow from the argument 
below for large m), 



a„ 



infr 



LceC",j|cj|=i ll^cjl 

where cTmin is the smallest singular value. For m large enough, there exist n points within 
- of n evenly spaced points u„. (Under the secondary hypothesis, this m clearly grows 
like 0{n'^).) Let Vg be the n x n Vandermonde matrix associated with these points, so 
that (under a certain ordering) 



Then 



\Vc\ 




> \\V,c\ 



We have 



K = K + -A, 



n 



where Vu is the Vandermonde matrix associated with u„ (i.e., a discrete Fourier transform) 



and ||A|| < 1. Thus ||V;c| 



which completes the proof. 



Vuc\\ + O(^). We know that 
ml VuC 

c6C",||c||=l 



\v- 



n 



D 



Lemma B.2. Suppose that, for all \x\ < 1 and e > 0, the smallest m such that min(|x — 
dm|) ^ e satisfies — = 0(e") for some a > 0. If f is analytic in the unit disk, then for 
m large enough the least squares approximation of f at the points d^ converges to f. 

Proof. Let 

Pn = {In, 0) 

denote the n x oo projection operator and let Em be the m x oo operator defined by 

Emf = /(dm,)- 
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Then we are approximating / by 
Furthermore, 

Pnf = V+EmPjPnf. 

We thus have the error 

f-f = f- PjPnf + P^V^E^P^PJ - /) 
= {I-PjV^Em){f-PjPnf). 

In other words, 

11/ — P^P„|| decays exponentially fast for any analytic /. The theorem follows since m 
grows at most algebraically with n. 

D 

Finally, we remark that the least squares system used in Algorithm |3] is not quite 
Vandermonde; it is actually of the form 

m 
k=l 

where / vanishes at —1. The logic of the preceding proofs still follow. To see this, define 
the n points Un as the points u„+i with the point —1 removed. Interpolating / at u^ 
by (z + 1, 2;^ — 1, . . . , z" — (—1)") will also interpolate / at —1, hence it is equivalent to 
interpolating / at u„+i. Thus the norm of the inverse of the relevant interpolation matrix 
at u„ is bounded by y/n + 1. 
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